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ABSTRACT 


For nearly two decades we have witnessed an intensive development 
Of a statistical methodology for assessing length of life and relia- 
bility of performance from empirical data. The initial stimulus for 
research on statistical problems in life testing and reliability came 
from the need to answer pressing practical questions which could not be 
treated by the existing statistical techniques. Because life and per- 
formance tests are so time consuming and expensive to run, it is a 
Practical necessity to terminate them as soon as possible. 

For the statistician this means developing estimation and decision 
procedure for data, which are severely curtailed in one way or another 
long before all items on test have actually failed. The 
estimation is more complicated when the data are truncated, i.e. when 
the observer loses track of some individuals before death occur. The 
product limit method of Kaplan and Meier is one way of estimating p(t) 
when the mechanism causing truncation is independent of the mechanism 
causing death. 

This paper proposes alternative estimators and compares them to 
the product limit method. A computer simulation is used to generate 
the times of death and truncation from a variety of assumed distribu- 
tions. No single estimator gives the best fit to the "true" distribu- 
tion of death under all situations. However, other estimators are 
shown to be better than the product limit estimator under all of the 


assumed situations. 
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I. INTRODUCTION 


Let the random variable T denote the time that elapses until an 
event occurs; the event may for example be an equipment failure, an 
individual's death, or the detection of a target. Denote by p(t) the 


probability of survival to time t, 
В = (Е) 


Picturesquely, T is called a lifetime, and p(t) is a survival pro- 


bÞability; 
ВЮ (Е) 15 the distribution function of T. 


In the medical field, one might wish to estimate the probability, 
p(t) that a patient survives t after a certain surgical procedure 
for cancer. In electronics, one wishes to estimate the probability of 
continuous failure-free operation of an equipment for time t. Іп the 
military, one might be interested in the probability of conducting a 
certain mission, under specified environmental conditions, without de- 
tection by the enemy. The event of interest may be a human death, 
equipment malfunctions, or sonar detection. Following Kaplan and Meier, 
Reference (1), this paper will refer to the event of interest as a 
"death". The test element in the sample may be a human, a radio, or 
a submarine. This paper will refer to the test elements as "individuals". 


Suppose that observed values of T are t a 


i 2 E SO that N 


lifetimes are observed. In this case an appropriate (unbiased) esti- 


= ОР Survival to time t is 





Under many circumstances complete lifetimes are not observed; censoring 
occurs at certains, х; , beyond which the life of an individual is not 
known. In such cases construction of an appropriate estimate of the 
survival probability is more difficult. In this paper various estimates 
of survival probability are studied when lifetimes are randomly censored. 
This means that censoring times are assumed to be realizations random 
variables independent of the actual lifetimes. 

The product-limit estimator of Kaplan and Meier, Reference (1), is 
an accepted method of dealing with the problem of censored data. This 
paper presents thirteen non-parametric estimators, including the product 
limit function. Censored data sets are simulated. The thirteen esti- 
mators are compared by examining their performance on the simulated data 


bases. 





те. THEORY 


There are two approaches to the empirical estimation of the survival 
probability, p(t): 

(1) one may use the observed fraction of survivors at arbitrarily 

selected times (step function estimator), or 

(2) one may focus attention on the times of the observed deaths 

(point estimator). 

The initial discussion is based on the assumption that all observa- 
tions are complete, i.e., it is assumed that all individuals remain under 
observation until their time of death. This initial assumption is for 
the purpose of simplifying the discussion. Then, later in this paper, 
the discussion is broadened to include incomplete data with observations 
of both death and censoring events. 

Survival Probabilities; No Censoring 


Let 0 = to ЕЕ -.. < t: < t < ... be a sequence of fixed 


T: 2 B 


EImes when 121€ T is a lifetime 
= > 
p(t.) pit ы.) 


and denote the conditional probability of survival to time tir given 


survival to tial by 
- >t. T> t. 
p(t; lti)? ae e, | il! 
iru) ШЕЛ) 
= LES = = Ex (1) 
pi iij! Pia 





ТЕ p(t,_,) = О , define DUI CHE - 0. 
Тһеп 
puo = РО ріс, |) = Se. ecole р Ney в. 
1 
= qt. 
IT pte, lt.) (2) 
nn 


where p = Р(Е; ) Ж БО” 1 . 


Observations on Uncensored Data at Fixed Times 

Let a sample of N individuals come under observation. They are 
all observed from birth (or the appropriate event defining time zero) 
until death. With the first approach, preselects a series of times, 
O < сі < t5 < ... before examining the observed time of death. In the 
medical follow-up example, one might select the times corresponding to 


exactly 1,2,3,... years after a surgical procedure for cancer. An esti- 


mate of the conditional probability of survival to eir given survival 


T EE 
Beeren) = NE (3) 


With N, elements were present at the beginning of the interval, i.e., 


at time t. 


1-1/ апа ri elements failed during the interval. 


For a set of data which is not censored, М; = М. - г. . Now 


replace probabilities by their estimates in (2): 


1 1 
лыд ~ М. Ss T 
Е “ri 
p(t.) = П ре 11) з 
j=1 ізі 7? 





a o ыс шиши 
т 1-2 НЕ 
і 
Ж; 
I J 
- 1 - 4 _ 


Now the estimate p(t.) is of the form 


at F eee E 
5 N (г, + к. + к.) 
p(t.) па, 


and this is the same as 


E S. 

p(t, ) ЕЕ 
where 5. is the number of the original group, of size N, that survive 
to ti. If it is assumed that the N individuals each have the survival 
probability p(t), and that they die independently, then Sir the random 
number that survive to time t: is binomially distributed, with 5. being 


a realized value of 5.. Then, considering the estimate as a random 


variable, 
~ 5, 
E 
and 
> N p(t.) 
E[p(t;)] т т: = p(t.) 
and 
В. ЕЦ Ер (е3) 
Var [p(t,) ] = — 


Consequently p(t.) is an unbiased and consistent estimate of pit.). 
This is true for every t, and can be shown to be true for all en 1=1, 


Zee ease will]. 


PO 





All of this indicates that the estimate suggested is likely to be a 
good one if the sample size, N, is large. 


Clearly p (t) € p(t, ,) . The survival probability, p(t), is thus 


1-1 
estimated at a fixed sequency of times. At each time point, t. being a 
1 
typical one, there are г, fewer survivors than at t. 1' where r. = 0,1, 
i- i 


2,...,М. Consequently a plot of p(t,) shows a non-decreasing step 


function, with downward steps of varying sizes at t.,t 


1 ТЗ Ф 


If the above times are close together, and if the time of death T, 
has a density function, then one can anticipate seeing values of r. that 
T 


are either zero or unity. 


The so-called second approach is really a limiting case of the first, 
as the time of intervals of measurement decrease indefinitely. Thus when 
a death (or loss) occurs it is only a single event. 

When no losses take place, the case now considered, the time t: ОЁ 
the ith death is a really a realization of a random variable, denoted by 
ti this means that p(t.) the probability of surviving t., is a random 


variable. It can be shown that the expected value of p(t.) is 





N-i+l : 
E [p(ti)] Б N+] , A М 
< < --- < = 
where E, 2, ES 


The derivation involves integrating 





Е . N! E ЕЕ N-i 
Elp(t,)] = Р(Е) * (1! [1-р (Е) ] Е ee 
О 
ee Neat) 
~ Nl 


by transformation from p(t) to x; see Cramér, Mathematical Methods of 


Statistics, H. Cramér, Princeton University Press, 1946. 


181 





Thus one is led to use 


№-1+1 
№+1 





p(ti) = (4) 


as an estimate of the value of p(ti), t: БЕЛЕ che ith timelof death: 
Expression (4) provides estimator of the survival function at times of 
observed deaths when there are no losses because of censoring. The 


estimator at the points E Е; < Е. < --- < Си can be connected by 


2 
straight lines, or a step function with step sizes 1/(N+1) may be used. 
The estimators of equation (4) give intuitively acceptable results. 
For example, if the sample consists of only a single individual (N=1), 
then death is equally likely to occur before or after the time at which 


the true (but unknown) survival function equals one half. Thus, the 


result of equation (4) is reasonable: 


т 1 
E[p(t,))] = 5 
The point estimates of the second approach always occur at the times of 
discontinuity forestimates from the first approach. For example, con- 
sider a data base (N=4) with deaths observed at times 1,3,4 and 7. The 


first approach gives the following step function estimate of the survival 


function: 
1572, Оне = 1 
IO қы 5-5 3 
S = o5 D ES d 
0/25 A. SIE << 
050 ER 7 
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The second approach gives the following point estimates 


p(o) » 1.0 
р(1) = 0.8 
р(3) = 0.6 
р(4) = 0.4 
р(7) = 0.2 


A graphic comparison of the results of the two approaches is given 


below: 





It is difficult to decide how to smooth out the step functions 
that result from the first approach. By connecting the tops of the 
"Stairsteps,"one places an upper bound on reasonable estimates. By 
connecting the bottom corners of the stairsteps, one places a lower 
bound on reasonable estimates. One might draw a smooth, decreasing 
curve that passes through all (or almost all) of the vertical faces 
of the step-function estimate. The second approach suggests method 


Of selecting a unique point on each of these vertical segments. 


13 





Incomplete observations 


When some of the observations are incomplete, equation 
(4) requires modification. The expected value of the survival function 
at the time of the first observed death may be written: 
М 


i 
Elp(t,)] = N +r (5) 


Here Ny is the effective size of the sample during the interval termi- 
nated by the time observed for the first death (0,ti). In the special 
Case of no censoring events, the value of Ni is unambiguous. It is 
equal to the initial sample size (ч, = ч). In this case equation (5) 
Yeduces to equation (4), 


Subsequent point estimates for t tare may be calculated 


2r 
iteratively: 


N. 





E(p(t,)] = * E[p(t, ,)] 


1 

№. +1 1 
1 

where to = 0 and N, is the effective sample size over the time interval 


о Thüs, 
д 


1-1! 
т E 
БР Е у (6) 
j=} J 


Variance of the estimators 

Kaplan and Meier, reference (1), give an expression for the exact 
calculation of the variance of step functions. They also discuss 
"Greenwood's formula," a large sample approximation that ignores terms 


Of order Им, ^. 
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Нега, reference (2), presence without derivation an expression for 


the variance of estimates using the second approach (point estimators): 


1 д 


н, ЇЧ. 2 
ce - IT c—— 





V(t,) = Var tE [p(t,)1) = ETT 


( 
ј=1 2-1 j 


The notation here follows that for the estimating equation (6). 


1:5 





TAE ESTIMATORS 


This section describes the nine non-parametric estimators and four 
jackknife estimators of the survival probability. It also describes 
the parametric estimator for an exponential decay function.  Exponential 
life distributions are the starting point for much of reliability theory 
and practice. The estimator derived from the exponential is regarded as 
"par" when the simulated data is based on an underlying exponential decay 
distribution for deaths. Thus, when deaths are exponentially distributed, 
the non-parametric estimators may be compared relative to each other, 
and they may be compared with the parametric estimator as a standard. 

A hypothetical data base, consisting of five individuals, is used 


to illustrate each of the estimators. This sample data base is as 


follows: 
Individual Time of Death Time of Truncation 
A Т = 
B Unknown (>2) 2 
C 3 = 
D Unknown (>6) 6 
E "i - 


The data have been arranged in time sequence of the death and trunca- 
tion events. In the medical example, the data might indicate that 
patients A, C and E were observed to die exactly 1, 3 and 7 years, re- 
spectively, after their surgery. However, B and D moved away or other- 


wise became unavailable to the observer at these times. Further, the 





cause of the unobservability is unrelated to the patient's health and 


life expectancy. 


A.  STEP-FUNCTION ESTIMATORS 
1. The First Estimator, IE 
p, (t) is a naive estimator; it is expected to perform poorly 
relative to the other estimators. P, only depends on the data from in- 
dividuals whose deaths are observed. It ignores any information from 
the partial lifetimes noted for the censored observations. p, (t) 15 


simply the fraction of individuals surviving to at least time t among 


those individuals whose time of death is known. It is a step function: 





А p, (t) 
E 1.0 
1-3 02667 
= 9333 
7-0 0.00 


The naive estimator, p, (t), takes no account of the successful 
survival intervals observed for the censored individuals. Therefore 
it is biased in a downward (pessimistic) direction. 

2. The Second Estimator, пы" 

p, (t) is the product-limit estimate. Kaplan and Meier, refer- 
ence (1), have shown that this is the maximum likelihood estimator. The 
Observed events, both deaths and truncations, are arranged in increasing 
order of occurrence: tir tore rtp where N is the number of individuals 
in the sample. 

Let p(ti) denote the cumulative probability of survival of an 
individual from time zero to time t;. Let p(t|t,) denote the conditional 


probability of surviving to time t(> ti), given that the individual has 


17 





already survived to time ti- Then, 
ous po (tus p(t.|t. ) (E-1) 
If we define E, = апа р(О) = 1, then 


T 
P,(t,) = Uy 9215455422 (E-2) 


The product limit estimator is in the form of equation (E-2) with 





N а 
== й If the event at t. is 
J truncation J 
B,te,lt,_,)= (E-3) 
N.-] 
2 If the event at t. is 
Ns a death 2 


Неке ш is the number of individuals observed surviving in the interval 
E Е < DE This formulation causes the product limit estimator to 
be insensitive to the exact time of the censoring events. 


The estimator is unity from time zero to the time of the first 


event, t reflecting the fact that all individuals in our example are 


I” 


observed to live until at least time ti: 


= If the event at time E, is a truncation, then the estimator 


remains at unity until at least time t Again, no deaths 


2- 

are observed in the sample before t5. 

= If the event at time Es is a death, then the estimator drops 
to (N-1)/N. This drop reflects the observed death of 1/N of 
the survival sample just prior to es 


Values of the estimator p, are calculated iteratively at successive 


values of Е.(і-1,2,... А). 
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The size of the survival sample declines as truncations and 
deaths remove individuals from observation. For the hypothetical data 


base listed above, one obtains: 





ч p, (t) 

0-1 то 
1-2 4/5 =0.8 
2-3 (4/5) x (3/3) = 0.8 
3-6 П С 0 555 
6-7 (VIE) So (VAL) = Bele 
7-со ey AL 1090 


The product-limit estimator explicitly accounts for the sur- 
vival of these individuals (up to the time of the last death before 
each censoring event). Thus, р. (Е) is a step function with a value 
that is not less than p, (t) for any value of t. If the sample contains 
no censoring, then p, (t) and p. (t) are identical. 

If the last event in the sample is a truncation rather than a 
death, then the modified data give the folloving estimate, i.e., 
individual E had disappeared from the observer at time 6.5 (so that 


the fact of E's death at time 7 is unknown). 


~ 





а p, (t) - Modified data 
0-1. то 

1-3 0.8 

3-6.5 (05538 


Since the time of the death for individual E is now unknown, 
one can only estimate that: 


98500) 4 Ы 595-5222 


T9 





If the analyst is willing to assume a functional form for the 
survival function, then he may calculate the manner in which the 
estimator DS) decreases to zero. However, the data alone are insuf- 
ficient when a strictly non-parametric ner is used. 

The product-limit estimator is a useful and intuitively appeal- 
ing method of dealing with incomplete observations. It has been wider 
used and studied. However, the product-limit has one disturbing 
characteristic: 

Most of the biological, physical or other causes of deaths pro- 
duce a survival probability that continuously decreases in time. 
It is, therefore, one may be a little uncomfortable estimating 
the survival probability with a step function. One is tempted 
to smooth the estimator to make it a monotonic decreasing func- 
ПОП ОБЕ. 

3. The Third Estimator, "p,(t)" 

p4(t) is a modification of p,(t). Like p, (t), it is a step 
function with discrete drops at those times corresponding to the observed 
deaths in the sample population. It may also be expressed as a product 


of conditional probabilities: 
i 

р = | р Е-4 

pj (t) a lc 00 (E-4) 


where the ty are the times of observed deaths and A is zero, The con- 


ditional probabilities on the right-hand side of Equation (E-4) differ 
somewhat from those in Equation (E-2): 


_ Nal 
= —— Е-5 
PE EE) 3 (Е-5) 
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Equation (E-5) differs from Equation (E-3) in the interpreta- 


tion of the numbers of individuals at risk. Here, the value of М. 15 


taken to be the average number of individuals observed surviving in 
the interval between the (k-1)st observed death and the kth observed 
death. The number of observed survivors decrease at intermediate times 


if events are censored, and hence the N, are not necessary integers. 


k 


The value of N. is regarded as the effective sample size for 


the interval from 8-1 БӘ Е In the sample data base shown above, 


Ke 
individual B is known to have survived from time 1 to time 2, or half 

of the interval between the first death at t=1 and the second death at 
t=3. Therefore, the estimator Рз treats individual B as half a parti- 
cipant in the interval between the death of individuals A and C. 


The effective sample size for this interval is then 3.5 


(n = 3 + = = 3,5) (full contributions from individuals C, D and E, 


plus a half contribution from B). For our hypothetical data base, the 





following values are calculated for Б. 








е Ра (6) 

col 5/5 = 1.0 
eS 4.5 x l0 = 0.8 
837 (2.5/3.5) х 0.8 = 0.571 
(7) (1.75/2.75) x 0.571 = 0.364 


The value of p, (t) can never be less than the corresponding value of 
p, (t). In the special case with no Censoring events the estimators 
p(t), p, (t) and p, (t) are identical. 

One might perturb the data by shifting the time of B's trunca- 
tion event down to 1+€ or up to 3-€, € arbitrarily small. The depend- 
ence of the estimator P4 upon the exact time of the censoring events 


may now be demonstrated. 
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Eor purposes of illustration, the time of the censoring event 


for individual B(t,) 1$ decreased from 2 to 1.1, then increased to 2.9. 





E p, (t£), 5-72 p(t), t, = 1.1 РАСЕ Ко 
0-1 120 1706 20 

Ш-3 0.80 0.80 0.80 

3-7 OZILI 02538 0.597 

(7) 0.364 0.342 0.380 


This example demonstrates an intuitively appealing characteris- 
tic of the estimator, P,- As the total observed survival time increases 
for the individuals in our sample (with deaths held constant), the value 
of the estimating function increases over at least a portion of its 
range. 

We may safely assume that the true survival function eventually 
tends to zero with time, since no physical or biological system lives 
forever. However, there are no observations on the survival of indi- 
viduals beyond time 7. The data only indicate that our step-function 
estimator drops to a value of .364 at t=7, but the nonparametric esti- 
mator gives no information about the survival function's subsequent 
decline from .364 to zero. However, the data alone are insufficient 


when a strictly nonparametric estimator is used. 


B. POINT ESTIMATOR 

As mentioned above, the estimators Pye p, and De are somewhat unde- 
sirable because they give step-function estimates for a continuous 
Survival function. The next three estimators D Р. апа Бе are modifi- 
cation of the first three. Again they provide estimates of the survival 


function only at those points in time that corresponde to observed deaths. 
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These estimators are specified by Equations (E-2) and (E-4), except for 
a substitution of the term (N+1) in place of (N). 

Since the point estimators have rigorous definitions at only dis- 
crete points in time, it is necessary to offer an interpolation rule. 
That is, we need a method of "connecting the dots." The method proposed 
here is to assume that the survival function declines in a piece-wise 
exponential decay between the discrete points in time. This procedure 
is equivalent to assuming that the hazard function is essentially con- 
stant between a consecutive pair of the discrete times, but that the 
hazard varies from one time period to the next. Such an assumption is 
intuitively acceptable unless one suspects violent fluctuations in the 
Mazard function. 

Пе Еве тако "А (Е)" 

p, (t) is analogous to p, (t) in that only those individuals ob- 
served to die are included in the sample. These two estimators are naive 
because they suppress all data from the survival times of individuals 
terminated from observation by censoring. 

These estimates, i.e., p} (t) апа Б, (©), tend to ignore informa- 
tion from the more long-lived individuals in the sample, and they may 
be expected to give biased estimates of the survival function. 

The point estimator p,(t) gives the following values with sample 


data base presented earlier in this section. 


Interpolation 
е Dre e S 
Qt’ &n(3/4) 
2 3 EL „да (2/3) 
(5) x 0.75 (De 
O) 
E 9.5 ШЕ 6 
(5 хо. 5 | 
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The interpolation for connecting the dots are as follows: 


interpo- 
lation 





2. The Estimator, AE E 
The estimator P, (t) similarly corresponds to the product-limit 
estimator p,(t). These two estimators use information from the indi- 
viduals on whom there are censored observations. Ра, like 55, does 
not exploit information about that portion of the censored observation 


after the death event (o£ some other individual) preceding the censor- 


ing event. 
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For our hypothetical data base the following values are 


calculated for PEDE 


Interpolation 

t Ps (t) t Ps (t) 
O 1.0 
1 MSS 0-1 ео 

3 5 = in CD 
3 (2) х (ва = 0.625 1-3 (=) е 

n 5 =з га G) 
7 (5)х 02625 = 0.312 3-7 (а) е 


Whenever censored observations are present, the estimator E never 
exceeds Bs (t). 

Rox pa tt), the value of N, is taken to be the number of surviv- 
Е аце 15 in the sample just before the observation of the ith 
death. This value is smaller than the number of surviving individuals 
ШЕЕ селе (1-1)st death if any truncation events occur in the 
interval. In fact, N, is the smallest number of surviving individuals 


observed at any time during the interval (t. ti). Thus р. mignt be 


als 
expected to introduce a bias by using values of N, that are, on the 
average, too small. However, this bias would be much less severe than 
the bias anticipated for the estimator p,(t). 

The estimators p, and Ba are insensitive to the precise times 
of the Censoring events. A change in the time of the censoring event 


for individual В бо 1+Е to 3-€, € arbitrarily small, does not alter the 


estimates from ЕЙ апа РЕ given in the preceding paragraph. 
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rhe Estimator, EAS 
The estimator pc (t) corresponds to Ба (6) by accounting for all 
of the survival time for the truncated observations. For our hypo- 


thetical data base, the following values are calculated for pg (t): 


Interpolation 
t pg (t) t pe (t) 
0 WO ol at’ &n(5/6) 
1 5/6 = 0.833 t-1 Sa 
=== {п (=== 
Б БО но 
E LI х = 0, 2 = 
(4-5) X 0.833 0.648 l-3 (e 
t-3 1.%5 
1.75 4 x (5-75) 
7 (5-75) ЖОО б48 = 0,412 3-7 (0.648)е 


The estimator p, (t) is based on the average number of surviving 
individuals noted in the various time intervals. These estimators give 
part credit for individuals whose lifetime is censored in mid-interval. 
The value of М. for pc (t) is an unweighted time average. If the obser- 
vation of an individual is truncated after 23$ of the interval has 
elapsed, then that individual contributes a value of 0.23 to Ni. Indivi- 
duals who are observed to survive the entire interval, and the individual 
whose death terminates the interval each contribute a value of 1.0 to N,- 
This interpretation of the effective sample size is approximate if the 
hazard is approximately constant over the interval. If the hazard function 
changes markedly within a time interval containing censored events, then 
this interpretation of the effective sample size is biased. Therefore, 
the procedure of determining the value of N, for the estimator Bao) is 


based on the implicit assumption that the survival function is locally 
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exponential. If the hazard function may be assumed to vary slowly over 
each of the time intervals (t. ur ti) then Р, would appear to be biased 
on an acceptable approximation. 

The estimator Pe: like Pr depends on the precise times of all 


deaths and censoring events. 








" pg (CO), t, =2 Pe (t), t, = 1.1 Ре, Е, = 2.9 
0 1.0 1.0 150 
1 5/6 = 0.833 5/6 = 0.833 5/6 = 0.833 
35 3005 E 3.95 _ 
3 (7) x 0.833 = 0.648. Caos) x 0.833 = 0.628. (7-95) х 0.833 - 0.665 
ШЕЛЕ _ 7/5 B ШЕГЕ B 
7 (6255) х 026498 = 0.412. 537) x 0.628 = 0.399. (555) х 0.665 = 0.423 


This illustrates that an increase (or decrease) in the total observed 
survival time causes an increase (or decrease) in the estimate pe over 
at least some of its time range. 

If the last event is a censored, and not an observed, death, 
these estimators also require definition for the time period starting 
with the time of the last death and ending with the time of the final 
censoring event. 

The method proposed here for p, (t) and p, (t) is to continue the 
exponential function used in the interval terminated by the time of the 
last death. This procedure can be illustrated with the modified data 


base used above in the discussion of P. and - 


C. THE BAYESIAN ESTIMATORS 
Consideration is next given to quasi-Bayesian estimators based on 
a uniform prior distribution on the unit interval. Let аа Бе the 


true survival times of N individuals which are censored on the right by 


N follow-up times Чуге Хы. It is assumed that the X; are independent, 


27 








identically distributed random variables with common distribution p(t) 


and we wish to estimate the survival function 
НЕ) ET - t) 


However, we only have available the data, 


СУ 


min do A Y.) 
i i 


I Y, 
х= al 


Qu ЈЕ X; 2 Yi, Шея 


H 
En 

О» 
| 


O, then 2. is called "a loss", and if 


О» 
| 


= 1, then 2. is called "a death". 
Then р,19, = oj) = р, 1х; zoom] eue) abs 


The maximum likelihood estimator for p(t) is 


p(t) == where 5= # 6 


is the number of successful tests, s has the binomial distribution. 


T9 


N S ‚U 
P (S|p) `2) DES ЕЕ, ОЕ 


eae) 


Е 
m 


The joint density of s and p is 


Ш--5 
р 2; С) р°(1-р) О р < 1, 5-0,1,... М. 


The marginal for s 15 


1 
1 5 Пт 5 5! (М-5)! 1 


О 
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for s=0,1,...N. Thus, averaging over the values of p, all of which are 


assumed to be equally likely, the values of s are equally likely to occur. 


The posterior for p then is 


l (n2) 


ВР (5162) —— а. _М-е 
жие СТ (3+1)Г (М-$+1) р (1-р) Oele 


а beta density with parameters s+1 and N-s+1. The mean of the posterior 
ig (s+1}| (N+2) ала the modal (maximum value) of the posterior is s/N; thus 


the Bayes estimate of p (given s survivers occur in the sample of N) is 


Then, equation (C-1) yields a step function and also has shown that the 
uniform prior has the effect of adding two individuals to the popula- 
tion at risk with one dying at time zero and the other essentially 
inmortal. 

The Bayesian estimators based on a uniform prior distribution on 
the unit interval are denoted 0100) Е BU and B 4t) , that correspond, 
respectively, to the estimators p,(t), p,(t) and p4(t). The sample 


data base thus gives the following estimates of the survival function: 





t Pt p,5(€) p34(€) 

0-1 4/5 = 0.8 6/7 = 0.857 6/7 = 0.857 

ПУБ = ос > x 0.857 = 0.714 (2) x 0.857 = 0.714 
5 

3-7 2/5 = 0.4 (9) x 0.714 = 0.536 2-5, ТИЕ 

D ae (=) х 0.536 = 0.268 CLA x 0.556 = 0.354 
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At the time of the final event (whether a death or a truncation), 
these step-function estimators drop to some positive value. Again, 
we have no data to indicate how the survival function proceeds to zero 


at subsequent times. 


D. THE JACKKNIFE ESTIMATOR 
We will assume that we observed, or have generated in a simulation, 
a survival probagility ЫЕ, j=1,...,n, from various sample sizes. 
Furthermore we have some parameter or characteristic E of the 
sample size which we wish to estimate with an estimator P(t). The 
jackknife estimator p(t,n) described below is an approximately unbiased 
estimator of ENE A modification of it has other useful properties. 
p_, (t,n-1) is the estimator from the sample of n of the X;'s with the 


ith value deleted from the sample. 


Se a) (п-1) Б. (с,п-1) = eM 


p; (tn) AD СЕТА =. 


1 T 


~ Й: ~ 
Еј ји = = EE) 


MEEI TS 
i es 


ПЕ 1 
Еће DEC called the PSEUDO-values. 

The PSEUDO-values can be used to obtain variance estimates of D(t,n) 
and to set approximate confidence limits, using Student's t. 

The idea is that the PSEUDO-values will be approximately indepen- 
dently and normally distributed. The jackknife estimator P(t,n) тама 
sample average so we form an estimate 52 of its variance given by 


ptn) 
the following relationship (Miller, 1974): 


P = 2 
Бр. (кип) - ZB, (tun) 





n-1 


N 
|o 
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This procedure is particularly useful if the number of data points 
is small, but it must be used with care. Note, that the estimator Е 
is designed to eliminate a = bias term in the estimator p(t,n). Of 
course the computational aspects of the complete jackknife can be quite 
onerous, especially if p(n) were, say, a complicated maximum likelihood 
estimator. Miller, reference (4) has shown that the product limit 
estimator is its own jackknife. 
Logistic Transformation 

Although one can legitimately jackknife the Kaplan-Meier estimate 
directly, there is some reason to believe that a preliminary transforma- 
tion will give improved results. Consequently, consider the transforma- 


Elon 


and notice that where the range of p(t) is from zero to unity, the above 
transformation makes the range of 2 run from -™ to о. The procedure 
utilized will be as follows. 
(A) Compute the overall estimate at a time point t, using all N data 
points, and using a "continuity" correction that has the effect 
of removing the effect of a zero in the logarithm (see D.R. COX, 


Analysis of Binary Data, Methuen Monograph): 


D (t) + a 
e A N 2M 
ey (0) + y 


(B) Compute the 2-values by leaving out each data point in turn when 


eom БЕЈ: for i-1,2,...,N. 


~ 1 
Е) + 57ү 
Q E N-1,-i 2 ) 


N-1,1 ~ ib 
— + €———— 
E E 
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(С) Form the pseudo-values 


(О) Compute 2, 5 


(E) Put approximate confidence (1-a)°100% limits on E[2] as follows 


L< E(%] <H 
=> SS 
2 +(-) t (N=1) Y a 


where Н (1) 1-0 


(F) Transform bash to obtain 


, and 
1+е 1+е 


The true value, p(t), should be enclosed between these levels for 
roughly (1-a) 100% of all samples. The coverage properties of this pro- 
cedure will now be checked by simulation: succesSive samples of size N 
will be selected, the jackknife limits H and L will be computed for each, 
¿E H 

< < 

z Á P(t) < 


l+e l+e 
illustrating performance are given subsequently. 








and a check will be made as to whether or not. Tables 
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D E is the logistic transformation estimator from the sample n 
, 


of the Х.'5 with the ith value deleted from the sample. 


z ЈЕ 


N-1,-i = 1 
ES Е 
FN-1,-i 2 (1-1) 


3, 
Di, 


= Ол -3.04|7-3.04| -1.89 





N 
|| 
e 
Mo 
| 
e 
¡A 
> 


= N-1,-i 
р EU = NE 
- Me m, - (N-1) p un 
l-pu (t)* 2N l-Py-1,-i © t 2 (8-1) 


z, m) are called PSEUDO-values of logistic transformation, the 
following values are calculated: 


2. : 
L 





i 


l 2 3 4 5 

t 

t, то 158 | > 198 | 2 198 
t, 2.198 | 2.198 | 2.198 | 2.198 
t. omeec |-3.3]4 | 2.446 | 29446 
e 0.606 |-3.314 | 2.446 | 2.446 
t. -3.0626|-3.0626|-3.0626|-3.0526|-7.162 
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Average of the pseudo-values 


Ni 
|| 
“ік 
| Ma 
N 


1 


Invert to find jackknife estimator of logistic transformation 


1 
m — + = 
ES 
uz ee 
С alo” ш 
l + e $ 


called the jackknife estimator 
of logistic transformation 


Variance of the 2; 


2 
5 ‚ = Аш ш = = 


The following values are calculated: 







ti 0.5484 0.646 
0.5484 0.646 
0.0568 5106 
0.0568 0:576 


The jackknife estimator for estimating variability and giving confidence 







interval. 
Tukey, reference (3) has suggested that in the jackknife procedure 
we consider the pseudo values 2, (N) as approximately independent and 


identically districuted and consequently, since Z is an average of 
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the ds (N), proceed as if 


1, 


о 


М 
1 124% 
Ga Zen) 


i m3 


1=1 
has t-distribution with N-1 d.F. 
If the 2. are approximately normal variates (Miller has shown) 


confidence bands for the unknown p(t) are given, as for the mean of 


any normal variate when estimated from sample size n. 


© 
= 2 
2 + = ti-a/2 (N-1) (р- 1) 
ше. 1 
ЕЕС Р(Е) + ZN, _ 5, 
z - — t ШЕ п авас =) СЕ. (N-1) 
im COME 1-5 (8) + э в E 
Tn) =Z = A *i-a/2 
cin) = 2+ E 1-0/2 
1. L(N) Í 1. L(N) i 
(1+ —)е— el ea (1+ =e 
2N 2N < p(t) < 2N Е AN 
1 + FENN) l + CLR 


The following values are calculated: 


517572 ='2. 776 
Lower Int. Upper Int. 
10 

10 

1.0 


ШЕШ? 


ti tj Af e] ст ја 
сл fa jw [M IH 


0.14 
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The basis for this leap of the imagination seems to be that ЈЕ 
2 = Х = x then the procedure for obtaining confidence intervals 
using equation (D-1) and pseudo-values is the same as the procedure 


using jackknife. Then if Xy = 2 апа 


п 
= l 
2 = = > Z. we have 
ne i 
1=1 
Ze = NY = (n-1 2 
- М N-1,-i 
М 
{ = х.) RES 
s 
=NX - (N-1 
B N N -1 
N N 
о р ХЕ. Хх 
j=1 а - 
Thus the pseudo value 
== 18 п — 
= and = о, 
1 і no a n 
1=1 
The pseudo values are independent if 2 =. and they are normal if 
x. is normal. 
E.  PARAMETRIC ESTIMATOR, IDEM 
This paper considers one additional estimator, denoted p, (t). e 


is a parametric estimator. Therefore, it is not really a competitor to 
the thirteen non-parametric estimators considered here. In general, a 
parametric estimator would not be used if the functional form were re- 


garded as unknown. Similarly, a non-parametric estimator would not 
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normally be used if the survival function were strongly suspected to 
have a specified form. 
р. (©) is the well known maximum likelihood estimator for the expo- 


nential distribution: 


x A 
P- (t) = е 


b t 
T = ———————————— 
where number of observed death 


In our sample data base, the total observed survival time is 19, and 


three deaths are observed. Thus, 


; БК, =1+2+3+6+7= 19 


i 
E- 3 
B 
= -3t/19 
and Р- (Е) = е / 


Calculations for selected times of interest yield the following esti- 


mates: 


р- 0) - 1.0 

p- (1) = 0.854 
р- Ө) - 0.623 
к - 02331 


The thirteen non-parametric estimators are compared for a variety of 
generating distributions for both the death mechanism and censoring 


mechanism. 
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ТУ. INSTRUCTIONS FOR USING PROGRAM 


INPUT 

Each input card bears nine variables. The distribution of time of 
death is entered in the first set of (five) columns, the censoring dis- 
tribution is entered in the second set of (ten) columns, a parameter of 
the censoring distribution is entered in the third set of (ten) columns, 
the number of replication is entered in the fourth set of (five) columns, 
the number of the event is entered in the fifth set of (five) columns. 
For the purpose of all print output used code "0" and "1" in the sixth 
set of (five) columns, the seed number is entered in the seventh set of 
(five) columns, after the card giving the time of the last event of a 
data set, a card with "0" or "1" in the column 50 is inserted, i.e., the 
"О" indicating more data sets to follow and "1" indicating the last data 
sets and t value is entered in the ninth set of (eight) columns. 


The distribution of timeof death and of censoring time used code as 


follows: 
Code Type of Distribution 
1 Uniform 
2 Exponential 
3 Delta function 


OUTPUT 
The output lists: 
1) the time of each observed failure 
2) estimated survival probability at that time 


3) the variance of that estimator 
4) result of goodness fit 
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а) 


Ы) 


с) 


mean error 
mean absolute error (ABS) 


root-mean-square error (RMS) 


5) total number of observed death 


6) confidence interval at particular time 


Definition of Fortran Variables 


NDIE 


NTRUNC 


XTRUNC 


NREPL 


NEVENT 


NWRITE 


NEND 


the distribution of time of death 

the distribution of censoring time 

the parameter of the distribution of censoring time 
number of replication 

number of event 

write all output or partial output of simulation 
indicate more data sets or last data set 

t statistic value 

the estimator, p, (t) 

the estimator, p, (t) 

the estimator, p, (t) 

the estimator, p, (t) 

the estimator, Ps (t) 

the estimator, p, (t) 

parametric estimator 25:202 

jackknife estimator of logistic transformation of p,(t) 
jackknife estimator of logistic transformation of p, (t) 
jackknife estimator of logistic transformation of Рё (Е) 
Bayesian estimator of p, (t) 

Bayesian estimator of p, (t) 


Bayesian estimator of p, (t) 
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P14 
(т.т) 


SBAR 


Var 


jackknife estimator of logistic transformation of p, (t) 
PSEUDO-value 

average of pseudo-value 

variance of estimator, p(t) 

variance of jackknife estimator 

mean of goodness fit 

absolute mean of goodness fit 

root mean square error 
upper confidence interval of 
lower confidence interval of 
upper confidence interval of pg (€) 
lower confidence interval of pg (€) 
upper confidence interval of pg (€) 
lower confidence interval of pg (t) 
upper confidence interval of 


lower confidence interval of 
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„ (Е)) and jackknife estimator of 
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To compare RMS with product limit (p 
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Computer output of the fourteen estimators 
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У. RESULTS OF THE SIMULATION 


This paragraph presents graphical comparisons of the 13 estimators 
based on simulated data. Each comparison is based on 1000 replications 
of a simulated data base. The bias and RMS error (square-root of mean- 
squared error) of each estimator depends on the parameters that control 
the simulated data base. No single estimator dominates all others 
under all conditions. 

The bias and RMS errors of the estimators depend on several factors: 
(A) The sample size (NEVENT) of individuals under observation at time 
zero affects the accuracy of the estimators. In general, a larger sample 
Size leads to a better estimate than a smaller sample. Values of NEVENT 
Selected for simulation are 5, 10, 25, and 50 (plus one simulation with 
NEVENT = 100). 

(B) The distribution of times at which the observations are censored 
(unless the individual dies earlier) affects the performance of the 
various estimators. This distribution is particularly important in con- 
junction with the distribution of lifetimes (do most individuals die 
before censoring is likely?, are deaths and censoring events about 
equally likely at all times?, are most observations censored before death?). 
Three types of distributions are assumed to underlie the censoring mech- 
anism: 
(1) Some of the samples are generated on the assumption that no censor- 
а оесико. 
(2) Some samples are generated from a uniform distribution of times of 


censoring. 
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(3) Other data bases are generated from an exponential distribution of 
censoring times. 

(С) The distribution of lifetimes (ignoring the possibility of censor- 

ing) also affects the performance of the various estimators. Two types 

of distributions are assumed to underlie the death mechanism: 

(1) Some of the samples are generated from a uniform distribution of 
lifetimes. 

(2) Other data bases are generated from an exponential distribution of 
lifetimes. 

If a uniform distribution of lifetimes is selected, its range is 
always over the interval from time O to time 1. If an exponential dis- 
tribution is selected, it always has a mean lifetime of l. The distri- 
butions of truncation times (uniform or exponential) have parameters 
E5607, ./5, l, 1.333, 1.5, 2 and 4. A wide"variety of samples 
may be simulated by mxing various pairs of distributions (for censoring 
times and deaths). Since the time units are arbitrary, the restriction 
on mean lifetimes is irrelevant. 

The true value of the survival function is, p(t), and the form of 
this function affects the relative performance of the 13 nonparametric 
estimators. For example, the Bayesian estimator 215 (5) tends to be 
better as measured by square-root of mean-squared error than its counter- 


part (the product-limit estimator, p, (t)) for the time frame in which 
а 29 


However, the product limit estimator tends to be better for those times 


when p(t) is close to zero or unity. 
The point estimators, Ps (t) and pe (€) tend to be better than the 


product-limit estimator (p, (t)) for all time periods. The jackknife 
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estimators of logistic transformation (pg (t), Pot), (BE) РОБ а 


a 
estimators tends as same as its counterpart point estimators (p, (tO), 

Ps (t), pg (t)) for all time periods. And also the estimator formed by 
jackknifing the logistic transformation (p, 4 (€)) of the product limit 


estimator tends to be better than its counterpart product limit (p, (t)) 


for the time frame in which 
ИВ (ЕЛ 7 


However, the product limit estimator tends to be better for those times 
when p(t) is close to unity. Point estimators, Ps (t) and pg (t) tend to 


be same for the time frame in which 
От = p, (Е) < 0.9 


However, the Ps (t) tends to be better for those times when р (©) 15 

close to unity. The jackknife procedure may be validated, in an empiri- 
cal sense, by sampling experiments or computer simulation in the follow- 
ign manner. First, times of censoring and death are obtained by drawing 
random numbers from postulated distributions. Second, the jackknifed 
estimator of the logistic-transformed product-limit estimation is found, 
and confidence limits are computed by the method of Tukey, reference (3). 
Since the true value of survival function, p(t), is known, so is the 
theoretical value of A. The jackknife confidence intervals can be 

checked for coverage: if Ly Sa EH then the particular interval covers, 
while otherwise (if A < Ly er H < A) it does not cover. Finally, the 
above procedure can be repeated many times (say 1000) and the fraction of 
repetitions which contains the true value of A is recorded. This fraction 
of the coverage should desirably be close to (l-a), the nominal confidence 


level. 
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The jackknife confidence limit procedure can be said to be robust 
of validity, ref (7), if the actual coverage is close to the nominal 
coverage, 1-a, for a various distributions. Such seems to be true for 
large n (n > 50). However, the jackknife confidence limits do not cover 
accurately when the true value of p(t) is close to unity. 

The following tables illustrate confidence limits of jackknife 
method of product limit (р,(Е)). Many computer generated graphics are 


presented on the following pages to complete this section. 
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APPENDIX A 


ESTIMATORS FOR GROUPED DATA 


With grouped Censored data the definition of p(t/t. .) given 
by equation (5) does not hold unless the assumption is made that all 
truncations occur at the end of the time interval. If, on the other 
hand, it is assumed that all truncations occur at the beginning of 
At, the equivalent form of equation (5) is 


N. - а. - Y. 
1 





plt,/t,_,) > (G-1) 


ШШЕ а. 
т al 


with N, elements were presents at beginning of interval, i.e., at 


time t. 


1-1' Е: elments failed during the interval, and а, elements 


truncated fror the sample during the interval but prior to failing. 
As a hypothesis, assume that all aborts occur simultaneously somewhere 
within the time interval, so that r' failures occur prior to the 


truncations and time remaining r; - r' after the truncations. Then 


N, zer N. т а= г. 
E EN Sg p i e == шшш ААЬАН G-2 
pl d e H. HU TEL es gs 
l l L 1 


Thus, the value of p(t;/t; j) depends on when the truncations occur. 
It is assumed that this is not known for the grouped data case.  Never- 


) 


theless, it is possible to place limits on the value of p(ti/t. 4 


Since equation (G-2) always gives values between those of equation (5) 


and (G-1). Thus 


(G-3) 
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For average sample size approximation, a simpler expression from the 
point of view of computational ease may be derived by substituting 
a/2 for a in equation (G-1) giving 


ү = 5 - Y 


p(t.) = — < (G-4) 
- ge 5 


The equation (G-4) may be thought of as the result of assuming that 
the average number of elements in the time interval is the number at 
the beginning decreased by halí the number of truncations. 

Records are usually available to provide a fairly precise time the 
death events. Іп the medical example, the exact time of death is 
usually recorded in medical records required by law. Іп the equipment 
lifetesting example, the time of malfunction or failure is usually 
known very precisely if the results are catastrophic; and maintenance 
records give a reasonably precise time even if the failure is not 
critical to a larger system. Іп the military example, the event of 
interest is usually a sensor detection or some other action that is 
routinely recorded in a log book. 

Equaiton (G-4) is a modification to the product-limit estimator, 
Pos when the times of truncation are known only in grouped form. Herd, 
reference (2), suggests a similar modification to estimators using the 
second approach (р. ок ре) with aggregated truncation data. Illustrate 
results for this method based on the sample data base of the nain test 
are given below. Here, of course, we do not know that individual B 
dropped out of observation at time 2 and that individual D dropped out 
at time 6. e know only that the two truncations occurred in the 


imterval (1,3) and (3,7), respectively. 
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Product limit's modification is denoted by p,' (t) and Herd's modi- 


fication is denoted by p,'(t). 


Thelr results on the sample data base are as 
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